###############
# R-code for replication
# Beyond innumeracy: Measuring public misperceptions about immigration

## Author: Philipp Lutz
## Date: 13.09.2025

## load libraries
library(foreign)
library(readstata13)
library(ggplot2)
library(ggpubr)
library(scales)
library(ggmosaic)
library(coefplot)
library(hrbrthemes)
library(corrplot)
library(vioplot)
library(patchwork)
library(tidyr)
library(tidyverse)
library(dplyr)
library(MASS)
library(haven)

## Load data
mosaich <- read_sav("swissubase_2503_1_0/data/2503_MOSAiCH2023_Data_v1.0.sav")
# Ernst St�hli, Mich�le, Sapin, Marl�ne, Pollien, Alexandre, Ochsner, Michael, and Nisple, Karin. (2024). 
# MOSAiCH 2023. Measurement and Observation of Social Attitudes in Switzerland. 
# Study on National Identity and Citizenship and related topics. (Version 1.0.0)
# [Dataset]. Lausanne: FORS data service. https://doi.org/10.48573/ae9w-0w40
time_stamps <- read.dta13("mosaich_time_stamps.dta") # available upon request from FORS data service (Lausanne)
mosaich <- merge(mosaich, time_stamps, by = "IDNO", all.x = TRUE)

## Data preparation

# Prepare demographic variables
mosaich$gender <- mosaich$DEMO1v2

# Immigration attitudes (higher value = more anti-immigration)
mosaich$immigration_item1 <- as.numeric(mosaich$NIC7a) - 3
mosaich$immigration_item2 <- 7 - as.numeric(mosaich$NIC7b)  # good for economy
mosaich$immigration_item3 <- as.numeric(mosaich$NIC7c) - 3  # take away jobs
mosaich$immigration_item4 <- 7 - as.numeric(mosaich$NIC7d)  # improve society
mosaich$immigration_item5 <- as.numeric(mosaich$NIC7e) - 3  # native privileges

mosaich$immigration <- with(mosaich, immigration_item1 + immigration_item2 +
                              immigration_item3 + immigration_item4 +
                              immigration_item5)
mosaich$immigration <- ((mosaich$immigration * -1) + 20)/20 # scale to a range from 0 to 1

# Calculate accuracy of perceptions
true_values <- c(31, 5, 14, 66, 4, 35, 54, 4)
mosaich$NICS9a_wp2[mosaich$NICS9a_wp2>100] <- NA
mosaich$NICS11a_wp2[mosaich$NICS11a_wp2>100] <- NA

# Define column name vectors once
perception_columns <- paste0("NICS", 8:15, "a_wp2")
confidence_vars <- paste0("NICS", 8:15, "b_wp2")

# Single integrated loop for all calculations
for (i in 1:8) {
  perc_col <- perception_columns[i]
  conf_col <- confidence_vars[i]
  true_val <- true_values[i]
  
  # Check for NA once per iteration
  perc_na <- is.na(mosaich[[perc_col]])
  
  # Calculate all accuracy measures at once
  raw_diff <- mosaich[[perc_col]] - true_val
  inaccuracy <- abs(raw_diff)
  inaccuracy_dir <- raw_diff
  
  mosaich[[paste0("inaccuracy", i)]] <- inaccuracy
  mosaich[[paste0("inaccuracy", i)]][perc_na] <- NA
  
  mosaich[[paste0("inaccuracy", i, "dir")]] <- inaccuracy_dir
  mosaich[[paste0("inaccuracy", i, "dir")]][perc_na] <- NA
  
  # Process confidence variables
  numeric_conf <- as.numeric(sub(" .*", "", as.character(mosaich[[conf_col]])))
  conf_na <- is.na(numeric_conf)
  
  # Calculate confidence measures
  confidence <- numeric_conf / 10
  nonconfidence <- 1 - confidence
  adj_confidence <- log(confidence + 1)
  adj_nonconfidence <- log(nonconfidence + 1)
  
  mosaich[[paste0("confidence", i)]] <- confidence
  mosaich[[paste0("confidence", i)]][conf_na] <- NA
  
  mosaich[[paste0("nonconfidence", i)]] <- nonconfidence  
  mosaich[[paste0("nonconfidence", i)]][conf_na] <- NA
  
  # Log-transformed versions
  mosaich[[paste0("adj_confidence", i)]] <- adj_confidence
  mosaich[[paste0("adj_nonconfidence", i)]] <- adj_nonconfidence
  
  # Max error and normalized measures
  max_error <- ifelse(raw_diff < 0, true_val, 100 - true_val)
  mosaich[[paste0("max_error", i)]] <- max_error
  norm_inaccuracy <- inaccuracy / max_error
  
  # All guessing variables
  mosaich[[paste0("guessing", i)]] <- inaccuracy * nonconfidence
  mosaich[[paste0("dguessing", i)]] <- inaccuracy * nonconfidence  # duplicate
  mosaich[[paste0("guessing", i, "dir")]] <- inaccuracy_dir * nonconfidence
  mosaich[[paste0("guessing", i, "log")]] <- inaccuracy * adj_nonconfidence
  
  # All misperception variables (INTEGRATED - no second loop needed)
  mosaich[[paste0("misperception", i)]] <- inaccuracy * confidence
  mosaich[[paste0("dmisperception", i)]] <- inaccuracy * confidence  # duplicate
  mosaich[[paste0("misperception", i, "log")]] <- inaccuracy * adj_confidence
  mosaich[[paste0("misperception", i, "dir")]] <- inaccuracy_dir * confidence
  
  # All normalized measures
  mosaich[[paste0("scmisperception", i, "_norm")]] <- norm_inaccuracy * confidence
  mosaich[[paste0("scguess", i, "_norm")]] <- norm_inaccuracy * nonconfidence
  mosaich[[paste0("scinaccuracy", i, "_norm")]] <- norm_inaccuracy
  mosaich[[paste0("scinaccuracy", i, "_dir")]] <- inaccuracy_dir / max_error
}

# Calculate overall metrics using vectorized operations (no changes needed here)
# Overall inaccuracy (mean)
mosaich$inaccuracy <- rowMeans(mosaich[paste0("inaccuracy", 1:8)])

# Overall confidence (mean)
mosaich$confidence <- rowMeans(mosaich[paste0("confidence", 1:8)])

# Normalized measures
mosaich$inacc_norm <- rowMeans(mosaich[paste0("scinaccuracy", 1:8, "_norm")])
mosaich$inacc_norm_dir <- rowMeans(mosaich[paste0("scinaccuracy", 1:8, "_dir")])

# Sum of normalized inaccuracies
mosaich$scinaccuracy <- rowSums(mosaich[paste0("scinaccuracy", 1:8, "_norm")])

# Directional sum with specific signs (items 4 and 6 are subtracted)
signs <- c(1, 1, 1, -1, 1, -1, 1, 1)
dir_data <- mosaich[paste0("scinaccuracy", 1:8, "_dir")]
mosaich$scinaccuracydir <- rowSums(sweep(dir_data, 2, signs, "*"))

## and now the correct calculation that inaccuracy means an inaccuracy in a negative direction
# Create a temporary copy of the 8 relevant variables
temp_inacc <- mosaich[paste0("scinaccuracy", 1:8, "_dir")]

# Invert the values for variables 4 and 6
temp_inacc[[4]] <- -temp_inacc[[4]]  # 
temp_inacc[[6]] <- -temp_inacc[[6]]  # 

# Compute the row-wise mean
mosaich$inacc_norm_dir <- rowMeans(temp_inacc, na.rm = TRUE)

# prepare misperception and guessing variables
mosaich$misperception <- (mosaich$misperception1 + mosaich$misperception2 + mosaich$misperception3 + mosaich$misperception4 + mosaich$misperception5 + mosaich$misperception6 + mosaich$misperception7 + mosaich$misperception8)/8
mosaich$misperceptionlog <- mosaich$misperception1log + mosaich$misperception2log + mosaich$misperception3log + mosaich$misperception4log + mosaich$misperception5log + mosaich$misperception6log + mosaich$misperception7log + mosaich$misperception8log
mosaich$misperception_norm <- (mosaich$scmisperception1_norm + mosaich$scmisperception2_norm + mosaich$scmisperception3_norm + mosaich$scmisperception4_norm + mosaich$scmisperception5_norm + mosaich$scmisperception6_norm + mosaich$scmisperception7_norm + mosaich$scmisperception8_norm)/8

mosaich$guessing <- (mosaich$guessing1 + mosaich$guessing2 + mosaich$guessing3 + mosaich$guessing4 + mosaich$guessing5 + mosaich$guessing6 + mosaich$guessing7 + mosaich$guessing8)/8
mosaich$guessinglog <- mosaich$guessing1log + mosaich$guessing2log + mosaich$guessing3log + mosaich$guessing4log + mosaich$guessing5log + mosaich$guessing6log + mosaich$guessing7log + mosaich$guessing8log
mosaich$guessing_norm <- (mosaich$scguess1_norm + mosaich$scguess2_norm + mosaich$scguess3_norm + mosaich$scguess4_norm + mosaich$scguess5_norm + mosaich$scguess6_norm + mosaich$scguess7_norm + mosaich$scguess8_norm)/8

### calculate knowledge and good guessing
mosaich$knowledge1 <- (69 - mosaich$inaccuracy1)*(mosaich$confidence1)
mosaich$knowledge2 <- (95 - mosaich$inaccuracy2)*(mosaich$confidence2)
mosaich$knowledge3 <- (84 - mosaich$inaccuracy3)*(mosaich$confidence3)
mosaich$knowledge4 <- (66 - mosaich$inaccuracy4)*(mosaich$confidence4)
mosaich$knowledge5 <- (96 - mosaich$inaccuracy5)*(mosaich$confidence5)
mosaich$knowledge6 <- (65 - mosaich$inaccuracy6)*(mosaich$confidence6)
mosaich$knowledge7 <- (54 - mosaich$inaccuracy7)*(mosaich$confidence7)
mosaich$knowledge8 <- (96 - mosaich$inaccuracy8)*(mosaich$confidence8)

mosaich$guessinggood1 <- (69 - mosaich$inaccuracy1)*(1-mosaich$confidence1)
mosaich$guessinggood2 <- (95 - mosaich$inaccuracy2)*(1-mosaich$confidence2)
mosaich$guessinggood3 <- (84 - mosaich$inaccuracy3)*(1-mosaich$confidence3)
mosaich$guessinggood4 <- (66 - mosaich$inaccuracy4)*(1-mosaich$confidence4)
mosaich$guessinggood5 <- (96 - mosaich$inaccuracy5)*(1-mosaich$confidence5)
mosaich$guessinggood6 <- (65 - mosaich$inaccuracy6)*(1-mosaich$confidence6)
mosaich$guessinggood7 <- (54 - mosaich$inaccuracy7)*(1-mosaich$confidence7)
mosaich$guessinggood8 <- (96 - mosaich$inaccuracy8)*(1-mosaich$confidence8)


## Data diagnostics

# survey mode test (web vs. paper)

# T-test for accuracy
t.test(inaccuracy ~ mode_wp1, data = mosaich)

# T-test for confidence
t.test(confidence ~ mode_wp1, data = mosaich)

## Data analysis

# Figure 1: Histograms of immigration perceptions
true_values <- c(Immigrants = 31, Refugees = 5, Muslims = 14, Europeans = 66,
                 Unemployed = 4, Educated = 35, Criminals = 54, Undocumented = 4)
titles <- names(true_values)
vars <- paste0("NICS", 8:15, "a_wp2")

# Generate perception variables and plots
plots <- list()

for (i in seq_along(vars)) {
  var <- vars[i]
  title <- titles[i]
  truth <- true_values[i]
  
  # Convert and clean
  mosaich[[var]] <- as.numeric(mosaich[[var]])
  mosaich[[var]][mosaich[[var]] > 100] <- NA
  mosaich[[paste0("perception", i)]] <- mosaich[[var]]
  
  # Create plot
  p <- ggplot(mosaich, aes_string(x = var)) +
    geom_histogram(alpha = 0.5, binwidth = 4) +
    xlim(0, 100) +
    xlab("Share in %") +
    ggtitle(title) +
    theme_ipsum() + theme(plot.margin = margin(10,10,10,10)) +
    theme(
      panel.grid.major = element_line(size = 0.15, color = "grey85"),
      panel.grid.minor = element_line(size = 0.05, color = "grey90")
    ) +
    geom_vline(xintercept = mean(mosaich[[var]], na.rm = TRUE), color = "black", size = 0.5, alpha = 0.8) +
    geom_vline(xintercept = truth, color = "#C06000", linetype = "solid", size = 1.5) +
    geom_vline(xintercept = median(mosaich[[var]], na.rm = TRUE), linetype = "longdash", color = "black", size = 0.5, alpha = 0.8)
  
  plots[[i]] <- p
}

jpeg("figure1.jpeg", width = 12, height = 6, units = "in", res = 800)
ggarrange(plotlist = plots, ncol = 4, nrow = 2)
dev.off()

# Calculate mean perceptions for each variable
mean_perceptions <- sapply(vars, function(v) mean(mosaich[[v]], na.rm = TRUE))
# Compute mean deviation from the true value
mean_diff <- mean_perceptions - true_values
# Combine into a summary table
deviation_summary <- data.frame(
  Variable = titles,
  True_Value = true_values,
  Mean_Perception = mean_perceptions,
  Mean_Difference = mean_diff
)
print(deviation_summary)


## Figure A1: Histograms of inaccuracy of immigration perceptions

### Cumulative density plots

# Shared elements
category_labels <- c("Immigrants", "Refugees", "Muslims", "Europeans", "Unemployed", "Educated", "Criminals", "Undocumented")

color_palette <- adjustcolor(
  c("#005050", "#006030", "#C06000", "#B08000", "#003366", "#505A60", "#802000", "#D4B28C"),
  alpha.f = 0.6
)
names(color_palette) <- paste0("inaccuracy", 1:8)  # will be adapted in the function

# Plotting function
plot_cdf <- function(data, pattern, title, xlab_label, xlim_vals = NULL, ylim_vals = c(0, 1)) {
  cols_to_plot <- paste0(pattern, 1:8) # select variables
  
  long_data <- data %>%
    pivot_longer(cols = all_of(cols_to_plot), names_to = "Group", values_to = "Value") %>%
    mutate(Value = as.numeric(Value)) %>%         # Ensure numeric
    filter(!is.na(Value) & Value >= 0)           # Remove NAs and negative values
  
  # Adjust colors to match variable names
  colors <- color_palette
  names(colors) <- cols_to_plot

  ggplot(long_data, aes(x = Value, color = Group)) +
    stat_ecdf(geom = "step", size = 1.7, pad = FALSE) +
    {if (!is.null(xlim_vals)) coord_cartesian(xlim = xlim_vals) else NULL} +
    scale_color_manual(values = colors, labels = category_labels) +
    xlab(xlab_label) +
    ylab("Cumulative Density") +
    ggtitle(title) +
    theme_ipsum() +
    theme(
      legend.title = element_blank(),
      legend.position = c(1, 0.7),
      legend.justification = c(1, 1),
      legend.background = element_rect(fill = "white", color = "black", size = 0.5),
      panel.grid = element_blank()
    ) +
    guides(color = guide_legend(override.aes = list(shape = NA, linetype = 1, size = 1.7), ncol = 1)) +
    {if (!is.null(ylim_vals)) ylim(ylim_vals) else NULL}
}


# Generate plots
cdf1 <- plot_cdf(mosaich, "inaccuracy", title = "Inaccuracy", xlab_label = "Inaccuracy")
cdf2 <- plot_cdf(mosaich, "confidence", title = "Confidence", xlab_label = "Confidence")
cdf3 <- plot_cdf(mosaich, "dmisperception", title = "Misperception", xlab_label = "Misperception", xlim_vals = c(0, 20))
cdf4 <- plot_cdf(mosaich, "dguessing", title = "Guessing", xlab_label = "Guessing", xlim_vals = c(0, 50))

jpeg("figure2.jpeg", width = 14, height = 10, units = "in", res = 800)
ggarrange(cdf1, cdf2, cdf3, cdf4, ncol = 2, nrow = 2)
dev.off()

# Correlation between key variables
cor(mosaich$confidence, mosaich$inaccuracy, use = "pairwise.complete.obs")
cor(mosaich$confidence, mosaich$inaccuracy, method = "spearman", use = "pairwise.complete.obs")

cor(mosaich$misperception, mosaich$guessing, use="complete.obs") # 
cor(mosaich$misperception, mosaich$guessing, method = "spearman", use = "pairwise.complete.obs")

## Figure A1: Histograms of inaccuracy of immigration perceptions
create_plot <- function(data, inaccuracy_var, title, mean_val, median_val) {
  ggplot(data, aes_string(x = inaccuracy_var)) + 
    geom_histogram(aes(y = ..density..), color = "gray", fill = "gray", bins = 20) + 
    geom_density(alpha = .2, fill = "#005050", lwd = 1.2, color = "#005050") + 
    geom_vline(xintercept = mean_val, color = "#C06000", linetype = "dashed", lwd = 0.4) + 
    geom_vline(xintercept = median_val, color = "#C06000", linetype = "solid", lwd = 0.4) + 
    xlab("Inaccuracy") + ylab("Density") + 
    ggtitle(title) + theme(plot.margin = margin(10,10,10,10)) +
    theme_ipsum()
}

# Select the confidence variables
conf_vars <- paste0("confidence", 1:8)
# Calculate mean confidence for each item
mean_conf <- sapply(conf_vars, function(v) mean(mosaich[[v]], na.rm = TRUE))
# Put in a summary table
confidence_summary <- data.frame(
  Variable = conf_vars,
  Mean_Confidence = mean_conf
)
print(confidence_summary)


# List of accuracy variables and corresponding titles
inaccuracy_vars <- paste0("inaccuracy", 1:8)
titles <- c("Immigrants", "Refugees", "Muslims", "Europeans", "Unemployed", "Educated", "Criminals", "Undocumented")

# Calculate mean and median for each accuracy variable and store in a list
mean_vals <- sapply(inaccuracy_vars, function(var) mean(mosaich[[var]], na.rm = TRUE))
median_vals <- sapply(inaccuracy_vars, function(var) median(mosaich[[var]], na.rm = TRUE))

# Create a list of plots
plots <- mapply(create_plot, 
                data = list(mosaich), 
                inaccuracy_var = inaccuracy_vars, 
                title = titles, 
                mean_val = mean_vals, 
                median_val = median_vals, 
                SIMPLIFY = FALSE)

jpeg("figure_A1.jpeg", width = 12, height = 6, units = "in", res = 800)
ggarrange(plotlist = plots, ncol = 4, nrow = 2)
dev.off()

# Figure A2: Histograms of confidence in immigration perceptions
# Select the relevant columns
conf_cols <- paste0("confidence", 1:8)
# Compute means
means <- sapply(mosaich[conf_cols], mean, na.rm = TRUE)
# Compute medians
medians <- sapply(mosaich[conf_cols], median, na.rm = TRUE)

# create plots
titles <- c("Immigrants", "Refugees", "Muslims", "Europeans",
            "Unemployed", "Educated", "Criminals", "Undocumented")
plots <- lapply(1:8, function(i) {
  ggplot(mosaich, aes_string(x = paste0("confidence", i))) +
    geom_histogram(aes(y = ..density..), color = "gray", fill = "gray", bins = 11) +
    geom_density(alpha = .2, fill = "#005050", lwd = 1.2, color = "#005050") +
    geom_vline(xintercept = means[paste0("confidence", i)],  color = "#C06000", linetype = "dashed", lwd = 0.4) +
    geom_vline(xintercept = medians[paste0("confidence", i)], color = "#C06000", linetype = "solid", lwd = 0.4) +
    xlab("Confidence") + ylab("Density") +
    ggtitle(titles[i]) + theme(plot.margin = margin(10,10,10,10)) +
    theme_ipsum()
})

jpeg("figure_A2.jpeg", width = 12, height = 6, units = "in", res=800)
do.call(ggarrange, c(plots, ncol = 4, nrow = 2))
dev.off()

# Figure A3: Histograms of misperception scores
ipsum_colors <- c("#005050", "#C06000")  # Muted blue and orange

p1 <- ggplot(mosaich, aes(x=misperception1)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Immigrants") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception1, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception1, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p2 <- ggplot(mosaich, aes(x=misperception2)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Refugees") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception2, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception2, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p3 <- ggplot(mosaich, aes(x=misperception3)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Muslims") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception3, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception3, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p4 <- ggplot(mosaich, aes(x=misperception4)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Europeans") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception4, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception4, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p5 <- ggplot(mosaich, aes(x=misperception5)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Unemployed") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception5, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception5, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p6 <- ggplot(mosaich, aes(x=misperception6)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Educated") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception6, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception6, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p7 <- ggplot(mosaich, aes(x=misperception7)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Criminals") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception7, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception7, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

p8 <- ggplot(mosaich, aes(x=misperception8)) + 
  geom_histogram(aes(y=..density..), color="gray", fill="gray", bins=20) + 
  xlab("Misperception") + ylab("Density") +
  geom_density(alpha=.2, fill="#005050", lwd=1.2, color="#005050") + 
  ggtitle("Undocumented") + 
  theme_ipsum() +
  geom_vline(aes(xintercept=mean(misperception8, na.rm=TRUE)), 
             linetype="dashed", color="black", lwd=0.4) +
  geom_vline(aes(xintercept=median(misperception8, na.rm=TRUE)), 
             linetype="solid", color="black", lwd=0.4)

jpeg("figure_A3.jpeg", width = 12, height = 6, units = "in", res=800)
ggarrange(p1,p2,p3,p4,p5,p6,p7,p8, ncol=4,nrow=2)
dev.off()

# Figure A4: Histograms of guessing scores
# Calculate mean and median for each guessing variable
guess_cols <- paste0("guessing", 1:8)
means_g <- sapply(mosaich[guess_cols], mean, na.rm = TRUE)
medians_g <- sapply(mosaich[guess_cols], median, na.rm = TRUE)

# Titles for guessing variables
titles_g <- c("Immigrants", "Refugees", "Muslims", "Europeans",
              "Unemployed", "Educated", "Criminals", "Undocumented")

# Create all plots
plots_g <- lapply(1:8, function(i) {
  ggplot(mosaich, aes_string(x = paste0("guessing", i))) +
    geom_histogram(aes(y = ..density..), color = "gray", fill = "gray", bins = 20) +
    geom_density(alpha = .2, fill = "#005050", lwd = 1.2, color = "#005050") +
    geom_vline(xintercept = means_g[paste0("guessing", i)],  color = "black", linetype = "dashed", lwd = 0.4) +
    geom_vline(xintercept = medians_g[paste0("guessing", i)], color = "black", linetype = "solid", lwd = 0.4) +
    xlab("Guessing") + ylab("Density") +
    ggtitle(titles_g[i]) +
    theme_ipsum()
})

jpeg("figure_A4.jpeg", width = 12, height = 6, units = "in", res = 800)
do.call(ggarrange, c(plots_g, ncol = 4, nrow = 2))
dev.off()

## Figure A5: Density plots

# Define plot labels and color mappings
group_labels <- c(
  "Immigrants", "Refugees", "Muslims", "Europeans", 
  "Unemployed", "Educated", "Criminals", "Undocumented"
)

custom_colors <- c(
  "#005050", "#006030", "#C06000", "#B08000", 
  "#003366", "#505A60", "#802000", "#D4B28C"
)

# Helper function for plotting
plot_density <- function(data, prefix, x_label, title, xlim_vals = NULL) {
  long_data <- data %>%
    pivot_longer(cols = paste0(prefix, 1:8), names_to = "Group", values_to = "Value")
  
  if (prefix == "confidence") {
    long_data <- filter(long_data, Group != "confidence")
  }
  
  group_names <- paste0(prefix, 1:8)
  names(custom_colors) <- group_names
  names(group_labels) <- group_names
  
  p <- ggplot(long_data, aes(x = Value, color = Group)) +
    geom_density(alpha = 0.6, lwd = 1.7) +
    scale_color_manual(values = adjustcolor(custom_colors, alpha.f = 0.6), 
                       labels = group_labels) +
    xlab(x_label) +
    ylab("Density") +
    ggtitle(title) +
    theme_ipsum() +
    theme(
      legend.title = element_blank(),
      legend.position = c(1, 1),
      legend.justification = c(1, 1),
      legend.background = element_rect(fill = "white", color = "black", size = 0.5),
      panel.grid = element_blank()
    ) +
    guides(color = guide_legend(
      override.aes = list(shape = NA, linetype = 1, size = 1.7),
      ncol = 1
    ))
  
  if (!is.null(xlim_vals)) {
    p <- p + xlim(xlim_vals)
  }
  
  return(p)
}

# Generate plots
density1 <- plot_density(mosaich, "inaccuracy", "Inaccuracy", "Inaccuracy")
density2 <- plot_density(mosaich, "confidence", "Confidence", "Confidence")
density3 <- plot_density(mosaich, "dmisperception", "Misperception", "Misperception", xlim_vals = c(0, 20))
density4 <- plot_density(mosaich, "dguessing", "Guessing", "Guessing", xlim_vals = c(0, 50))

jpeg("figure_A5.jpeg", width = 14, height = 10, units = "in", res=800)
ggarrange(density1,density2,density3,density4,ncol=2,nrow=2)
dev.off()

## Figure A6: Cumulative density plots with normalised scores
group_labels <- c(
  "Immigrants", "Refugees", "Muslims", "Europeans",
  "Unemployed", "Educated", "Criminals", "Undocumented"
)

group_colors <- adjustcolor(
  c("#005050", "#006030", "#C06000", "#B08000",
    "#003366", "#505A60", "#802000", "#D4B28C"),
  alpha.f = 0.6
)
names(group_colors) <- paste0("prefix", 1:8)  # placeholder names

# Function to generate ECDF plots
generate_density_plot <- function(data, prefix, title, xlabel, xlim_vals = NULL) {
  # Match only the *_norm columns
  matched_cols <- grep(paste0("^", prefix, "[0-9]+_norm$"), names(data), value = TRUE)
  
  if (length(matched_cols) == 0) {
    stop(paste("No columns found for", prefix))
  }
  
  # Assign real names to colors
  names(group_colors) <- matched_cols
  
  long_data <- data %>%
    pivot_longer(cols = all_of(matched_cols), names_to = "Group", values_to = "Value")
  
  p <- ggplot(long_data %>% filter(!is.na(Value)), aes(x = Value, color = Group)) +
    stat_ecdf(geom = "step", size = 1.7, pad = FALSE) +
    scale_color_manual(
      values = group_colors,
      labels = group_labels
    ) +
    xlab(xlabel) +
    ylab("Density") +
    ggtitle(title) +
    theme_ipsum() +
    theme(
      legend.title = element_blank(),
      legend.position = c(1, 1),
      legend.justification = c(1, 1),
      legend.background = element_rect(fill = "white", color = "black", size = 0.5),
      panel.grid = element_blank()
    ) +
    guides(
      color = guide_legend(
        override.aes = list(shape = NA, linetype = 1, size = 1.7),
        ncol = 1
      )
    )
  
  if (!is.null(xlim_vals)) {
    p <- p + xlim(xlim_vals)
  }
  
  return(p)
}


# Now call the function
density1 <- generate_density_plot(mosaich, "scinaccuracy", "Inaccuracy", "Inaccuracy (normalised)")
density3 <- generate_density_plot(mosaich, "scmisperception", "Misperception", "Misperception (normalised)", c(0, 1))
density4 <- generate_density_plot(mosaich, "scguess", "Guessing", "Guessing (normalised)", c(0, 1))

jpeg("figure_A6.jpeg", width = 14, height = 10, units = "in", res=800)
ggarrange(density1,cdf2,density3,density4,ncol=2,nrow=2)
dev.off()


## Figure A7: Scatterplots of inaccuracy and confidence
mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy1) & !is.na(mosaich$confidence1))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy1, mosaich_nm$confidence1), mosaich_nm[,c("inaccuracy1", "confidence1")])
set.seed(123)  # Set seed for reproducibility
plot1 <- ggplot(mosaich_nm, aes(inaccuracy1, confidence1, color = misperception1, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE) + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Immigrants") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy2) & !is.na(mosaich$confidence2))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy2, mosaich_nm$confidence2), mosaich_nm[,c("inaccuracy2", "confidence2")])
set.seed(123)  # Set seed for reproducibility
plot2 <- ggplot(mosaich_nm, aes(inaccuracy2, confidence2, color = misperception2, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Refugees") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy3) & !is.na(mosaich$confidence3))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy3, mosaich_nm$confidence3), mosaich_nm[,c("inaccuracy3", "confidence3")])
set.seed(123)  # Set seed for reproducibility
plot3 <- ggplot(mosaich_nm, aes(inaccuracy3, confidence3, color = misperception3, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Muslims") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy4) & !is.na(mosaich$confidence4))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy4, mosaich_nm$confidence4), mosaich_nm[,c("inaccuracy4", "confidence4")])
set.seed(123)  # Set seed for reproducibility
plot4 <- ggplot(mosaich_nm, aes(inaccuracy4, confidence4, color = misperception4, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Europeans") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy5) & !is.na(mosaich$confidence5))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy5, mosaich_nm$confidence5), mosaich_nm[,c("inaccuracy5", "confidence5")])
set.seed(123)  # Set seed for reproducibility
plot5 <- ggplot(mosaich_nm, aes(inaccuracy5, confidence5, color = misperception5, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Unemployed") +
  scale_color_gradient(low ="#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy6) & !is.na(mosaich$confidence6))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy6, mosaich_nm$confidence6), mosaich_nm[,c("inaccuracy6", "confidence6")])
set.seed(123)  # Set seed for reproducibility
plot6 <- ggplot(mosaich_nm, aes(inaccuracy6, confidence6, color = misperception6, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Educated") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy7) & !is.na(mosaich$confidence7))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy7, mosaich_nm$confidence7), mosaich_nm[,c("inaccuracy7", "confidence7")])
set.seed(123)  # Set seed for reproducibility
plot7 <- ggplot(mosaich_nm, aes(inaccuracy7, confidence7, color = misperception7, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE)  + xlab("Inaccuracy") + ylab("Confidence") + 
  theme_ipsum() + ggtitle("Criminals") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

mosaich_nm <- subset(mosaich, !is.na(mosaich$inaccuracy8) & !is.na(mosaich$confidence8))
mosaich_nm$density1 <- fields::interp.surface(
  MASS::kde2d(mosaich_nm$inaccuracy8, mosaich_nm$confidence8), mosaich_nm[,c("inaccuracy8", "confidence8")])
set.seed(123)  # Set seed for reproducibility
plot8 <- ggplot(mosaich_nm, aes(inaccuracy8, confidence8, color = misperception8, alpha = 1/density1)) +
  geom_jitter(shape = 16, size = 5, show.legend = FALSE) +
  xlab("Inaccuracy") +
  ylab("Confidence") +
  theme_ipsum() +
  ggtitle("Undocumented") +
  scale_color_gradient(low = "#005050", high = "#C06000") +
  scale_alpha(range = c(.05, .25)) +
  stat_smooth(method = "loess", fullrange = TRUE, fill = "#C06000", color = "#C06000") +
  theme(legend.position = "none")

jpeg("figure_A7.jpeg", width = 12, height = 6, units = "in", res=800)
ggarrange(plot1,plot2,plot3,plot4,plot5,plot6,plot7,plot8, ncol=4,nrow=2)
dev.off()

## Figure A8: Scatterplot of inaccuracy and confidence (aggregated)
jpeg("figure_A8.jpeg", width = 10, height = 5, units = "in", res=800)
set.seed(123)  # Set seed for reproducibility
plot1 <- ggplot(mosaich, aes(x = inacc_norm, y = confidence)) + 
  geom_jitter(size = 3, color = "black", alpha=0.2) + theme_ipsum() +
  stat_smooth(method = "lm", fullrange = TRUE, fill="#C06000", color="#C06000") + xlab("Inaccuracy (absolute)") + ylab("Confidence") +
  geom_rug()
set.seed(123)  # Set seed for reproducibility
plot2 <- ggplot(mosaich, aes(x = inacc_norm_dir, y = confidence)) + 
  geom_jitter(size = 3, color = "black", alpha=0.2) + theme_ipsum() +
  stat_smooth(method = "lm", fullrange = TRUE, fill="#C06000", color="#C06000") + xlab("Inaccuracy (directional)") + ylab("Confidence") +
  geom_rug()
ggarrange(plot1,plot2)
dev.off()

## Figure 3: Average misperception and guessing by item
jpeg("figure3.jpeg", width = 6, height = 4.5, units = "in", res=800)
par(
  mar = c(5, 4, 0, 1),          # Margins (bottom, left, top, right)
  mgp = c(2, 0.5, 0),           # Axis title, label, and line positions
  tck = -0.01,                  # Tick mark length
  las = 1,                      # Horizontal axis labels
  font.main = 1,                # Plain font for title
  font.lab = 2,                 # Bold font for axis labels
  font.axis = 1,                # Plain font for axis text
  bty = "n",                    # Remove box around plot
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30",           # Axis label color
  fg = "gray30"                 # Foreground color
)

ipsum_colors <- c("#005050", "#C06000")  

# Calculate means for misperception and guessing
means <- c(
  mean(mosaich$misperception2, na.rm = TRUE), mean(mosaich$guessing2, na.rm = TRUE),
  mean(mosaich$misperception4, na.rm = TRUE), mean(mosaich$guessing4, na.rm = TRUE),
  mean(mosaich$misperception7, na.rm = TRUE), mean(mosaich$guessing7, na.rm = TRUE),
  mean(mosaich$misperception6, na.rm = TRUE), mean(mosaich$guessing6, na.rm = TRUE),
  mean(mosaich$misperception5, na.rm = TRUE), mean(mosaich$guessing5, na.rm = TRUE),
  mean(mosaich$misperception3, na.rm = TRUE), mean(mosaich$guessing3, na.rm = TRUE),
  mean(mosaich$misperception1, na.rm = TRUE), mean(mosaich$guessing1, na.rm = TRUE),
  mean(mosaich$misperception8, na.rm = TRUE), mean(mosaich$guessing8, na.rm = TRUE)
)

# Define category names
categories <- c("Refugees", "Europeans", "Criminals", "Educated", 
                "Unemployed", "Muslims", "Immigrants", "Undocumented")

# Create barplot without x-axis labels
bp <- barplot(
  means,
  names.arg = NULL,             # Remove default labels
  ylab = "Inaccuracy",
  border = NA,                  # No border around bars
  col = rep(ipsum_colors, 8),   # Use ipsum color palette
  beside = TRUE,                # Show bars side-by-side
  space = c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0),  # Adjust spaces between pairs
  ylim = c(0, max(means, na.rm = TRUE) * 1.1),  # Add padding to y-axis
  cex.axis = 0.8,               # Adjust size of axis labels
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30"            # Axis label color
)

abline(
  h = c(0, 5, 10, 15),      
  col = "gray90",               
  lty = 2                     
)

# Add custom x-axis labels (centered between pairs of bars)
text(
  x = colMeans(matrix(bp, nrow = 2)),  # Calculate midpoints between bars
  y = par("usr")[3] - 0.05 * (par("usr")[4] - par("usr")[3]),  # Position below bars
  labels = categories,          # Category names
  srt = 45,                     # Rotate labels for better readability
  adj = 1,                      # Adjust alignment
  xpd = TRUE,                   # Allow text outside plot area
  cex = 0.8,                    # Adjust text size
  col = "gray30"                # Match text color to theme
)

legend(
  x = max(bp) * 0.5,            # X-coordinate (right side of the plot)
  y = 18,                      # Y-coordinate (lower position)
  legend = c("Misperception", "Guessing"),  # Labels for the colors
  fill = ipsum_colors,          # Colors corresponding to the bars
  border = NA,                  # No border for the legend boxes
  bty = "n",                    # Remove legend box outline
  cex = 0.9,                    # Adjust text size in the legend
  text.col = "gray30",          # Match text color to theme
  ncol = 2                      # Arrange legend items vertically
)
dev.off()

## Figure A9: Average misperception and guessing based on normalized scores
jpeg("figureA9.jpeg", width = 6, height = 4.5, units = "in", res=800)
par(
  mar = c(5, 4, 0, 1),          # Margins (bottom, left, top, right)
  mgp = c(2, 0.5, 0),           # Axis title, label, and line positions
  tck = -0.01,                  # Tick mark length
  las = 1,                      # Horizontal axis labels
  font.main = 1,                # Plain font for title
  font.lab = 2,                 # Bold font for axis labels
  font.axis = 1,                # Plain font for axis text
  bty = "n",                    # Remove box around plot
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30",           # Axis label color
  fg = "gray30"                 # Foreground color
)
par(mfrow=c(1,1))

# Calculate means for misperception and guessing
means <- c(
  mean(mosaich$scmisperception4_norm, na.rm = TRUE), mean(mosaich$scguess4_norm, na.rm = TRUE),
  mean(mosaich$scmisperception6_norm, na.rm = TRUE), mean(mosaich$scguess6_norm, na.rm = TRUE),
  mean(mosaich$scmisperception7_norm, na.rm = TRUE), mean(mosaich$scguess7_norm, na.rm = TRUE),
  mean(mosaich$scmisperception2_norm, na.rm = TRUE), mean(mosaich$scguess2_norm, na.rm = TRUE),
  mean(mosaich$scmisperception3_norm, na.rm = TRUE), mean(mosaich$scguess3_norm, na.rm = TRUE),
  mean(mosaich$scmisperception1_norm, na.rm = TRUE), mean(mosaich$scguess1_norm, na.rm = TRUE),
  mean(mosaich$scmisperception5_norm, na.rm = TRUE), mean(mosaich$scguess5_norm, na.rm = TRUE),
  mean(mosaich$scmisperception8_norm, na.rm = TRUE), mean(mosaich$scguess8_norm, na.rm = TRUE)
)

# Define category names
categories <- c("Europeans", "Educated", "Criminals",  "Refugees", 
                  "Unemployed", "Muslims","Immigrants", "Undocumented")

# Create barplot without x-axis labels
bp <- barplot(
  means,
  names.arg = NULL,             # Remove default labels
  ylab = "Magnitude of inaccurate perception",
  border = NA,                  # No border around bars
  col = rep(ipsum_colors, 8),   # Use ipsum color palette
  beside = TRUE,                # Show bars side-by-side
  space = c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0),  # Adjust spaces between pairs
  ylim = c(0, max(means, na.rm = TRUE) * 1.1),  # Add padding to y-axis
  cex.axis = 0.8,               # Adjust size of axis labels
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30"            # Axis label color
)

abline(
  h = c(0, 0.05,0.1,0.15,0.2,0.25,0.3),       # 
  col = "gray90",               # 
  lty = 2                       # 
)

# Add custom x-axis labels (centered between pairs of bars)
text(
  x = colMeans(matrix(bp, nrow = 2)),  # Calculate midpoints between bars
  y = par("usr")[3] - 0.05 * (par("usr")[4] - par("usr")[3]),  # Position below bars
  labels = categories,          # Category names
  srt = 45,                     # Rotate labels for better readability
  adj = 1,                      # Adjust alignment
  xpd = TRUE,                   # Allow text outside plot area
  cex = 0.8,                    # Adjust text size
  col = "gray30"                # Match text color to theme
)

legend(
  x = max(bp) * 0.55,            # X-coordinate (right side of the plot)
  y = 0.29,                      # Y-coordinate (lower position)
  legend = c("Misperception", "Guessing"),  # Labels for the colors
  fill = ipsum_colors,          # Colors corresponding to the bars
  border = NA,                  # No border for the legend boxes
  bty = "n",                    # Remove legend box outline
  cex = 0.9,                    # Adjust text size in the legend
  text.col = "gray30",          # Match text color to theme
  ncol = 2                      # Arrange legend items vertically
)
dev.off()

# Figure A10
jpeg("figure_A10.jpeg", width = 6, height = 4.5, units = "in", res=800)
par(
  mar = c(5, 4, 0, 1),          # Margins (bottom, left, top, right)
  mgp = c(2, 0.5, 0),           # Axis title, label, and line positions
  tck = -0.01,                  # Tick mark length
  las = 1,                      # Horizontal axis labels
  font.main = 1,                # Plain font for title
  font.lab = 2,                 # Bold font for axis labels
  font.axis = 1,                # Plain font for axis text
  bty = "n",                    # Remove box around plot
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30",           # Axis label color
  fg = "gray30"                 # Foreground color
)

# Define a clean color palette (ipsum-style)
ipsum_colors <- c("#005050", "#C06000")  # Muted blue and orange

# Calculate means for misperception and guessing
means <- c(
  mean(mosaich$knowledge5, na.rm = TRUE), mean(mosaich$guessinggood5, na.rm = TRUE),
  mean(mosaich$knowledge1, na.rm = TRUE), mean(mosaich$guessinggood1, na.rm = TRUE),
  mean(mosaich$knowledge8, na.rm = TRUE), mean(mosaich$guessinggood8, na.rm = TRUE),
  mean(mosaich$knowledge3, na.rm = TRUE), mean(mosaich$guessinggood3, na.rm = TRUE),
  mean(mosaich$knowledge2, na.rm = TRUE), mean(mosaich$guessinggood2, na.rm = TRUE),
  mean(mosaich$knowledge4, na.rm = TRUE), mean(mosaich$guessinggood4, na.rm = TRUE),
  mean(mosaich$knowledge6, na.rm = TRUE), mean(mosaich$guessinggood6, na.rm = TRUE),
  mean(mosaich$knowledge7, na.rm = TRUE), mean(mosaich$guessinggood7, na.rm = TRUE)
)

# Define category names
categories <- c("Unemployed", "Immigrants", "Undocumented", 
                "Muslims", "Refugees", "Europeans", "Educated", "Criminals")

# Create barplot without x-axis labels
bp <- barplot(
  means,
  names.arg = NULL,             # Remove default labels
  ylab = "Inaccuracy",
  border = NA,                  # No border around bars
  col = rep(ipsum_colors, 8),   # Use ipsum color palette
  beside = TRUE,                # Show bars side-by-side
  space = c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0),  # Adjust spaces between pairs
  ylim = c(0, max(means, na.rm = TRUE) * 1.1),  # Add padding to y-axis
  cex.axis = 0.8,               # Adjust size of axis labels
  col.axis = "gray30",          # Axis text color
  col.lab = "gray30"            # Axis label color
)
abline(
  h = c(0, 10,20,30,40,50,60,70),       # 
  col = "gray90",               # 
  lty = 2                       #
)

# Add custom x-axis labels (centered between pairs of bars)
text(
  x = colMeans(matrix(bp, nrow = 2)),  # Calculate midpoints between bars
  y = par("usr")[3] - 0.05 * (par("usr")[4] - par("usr")[3]),  # Position below bars
  labels = categories,          # Category names
  srt = 45,                     # Rotate labels for better readability
  adj = 1,                      # Adjust alignment
  xpd = TRUE,                   # Allow text outside plot area
  cex = 0.8,                    # Adjust text size
  col = "gray30"                # Match text color to theme
)

legend(
  x = max(bp) * 0.4,            # X-coordinate (right side of the plot)
  y = 65,                      # Y-coordinate (lower position)
  legend = c("Knowledge", "Guessing"),  # Labels for the colors
  fill = ipsum_colors,          # Colors corresponding to the bars
  border = NA,                  # No border for the legend boxes
  bty = "n",                    # Remove legend box outline
  cex = 0.9,                    # Adjust text size in the legend
  text.col = "gray30",          # Match text color to theme
  ncol = 2                      # Arrange legend items vertically
)
dev.off()

mosaich$knowledge <- mosaich$knowledge1 + mosaich$knowledge2 + mosaich$knowledge3 + mosaich$knowledge4 + mosaich$knowledge5 + mosaich$knowledge6 + mosaich$knowledge7 + mosaich$knowledge8
mosaich$guessinggood <- mosaich$guessinggood1 + mosaich$guessinggood2 + mosaich$guessinggood3 + mosaich$guessinggood4 + mosaich$guessinggood5 + mosaich$guessinggood6 + mosaich$guessinggood7 + mosaich$guessinggood8
cor(mosaich$knowledge, mosaich$guessinggood, use="complete.obs") # 
# Calculate row-wise means for knowledge and guessinggood variables
mosaich$knowledge_mean <- rowMeans(mosaich[, paste0("knowledge", 1:8)], na.rm = TRUE)
mosaich$guessinggood_mean <- rowMeans(mosaich[, paste0("guessinggood", 1:8)], na.rm = TRUE)
# Calculate correlation between the two means
correlation <- cor(mosaich$knowledge_mean, mosaich$guessinggood_mean, use = "complete.obs")

# Print the result
print(correlation)

## correlation test with immigration attitudes
cor.test(mosaich$misperception, mosaich$immigration, method = "spearman") # 
cor.test(mosaich$guessing, mosaich$immigration, method = "spearman") # 

library(cocor)
# Clean data: remove NAs in any of the relevant columns
mosaich_clean <- na.omit(mosaich[, c("misperception", "guessing", "immigration")])
n <- nrow(mosaich_clean)
# Compute Spearman correlations directly
cor_mis_imm <- cor(mosaich_clean$misperception, mosaich_clean$immigration, method = "spearman")
cor_gue_imm <- cor(mosaich_clean$guessing, mosaich_clean$immigration, method = "spearman")
cor_mg     <- cor(mosaich_clean$misperception, mosaich_clean$guessing, method = "spearman")
# Test difference between dependent correlations with overlapping variable
cocor.dep.groups.overlap(
  r.jk = cor_mis_imm,  # 
  r.jh = cor_gue_imm,  # 
  r.kh = cor_mg,        # 
  n = n,
  alternative = "two.sided"
)

## decomposition analysis: inaccuracy and (non-)confidence
library("relaimpo")  # 
model_mis <- lm(misperception ~ inaccuracy + confidence, data = mosaich)
model_guess <- lm(guessing ~ inaccuracy + I(1 - confidence), data = mosaich)

# Compute relative importance
shapley_mis <- calc.relimp(model_mis, type = "lmg", rela = TRUE)
shapley_guess <- calc.relimp(model_guess, type = "lmg", rela = TRUE)

# Print results
print(shapley_mis)
print(shapley_guess)

# Visualizing results
mis_data <- data.frame(
  Variable = c("Inaccuracy", "Confidence"),
  Importance = shapley_mis$lmg
)
guess_data <- data.frame(
  Variable = c("Inaccuracy", "1 - Confidence"),
  Importance = shapley_guess$lmg
)

# Misperception Score Plot
mis_plot <- ggplot(mis_data, aes(x = Variable, y = Importance, fill = Variable)) +
  geom_bar(stat = "identity", show.legend = FALSE) +
  scale_fill_manual(values = c("#C06000", "#005050")) +
  labs(
    title = "Misperception score",
    y = "Relative importance",
    x = NULL
  ) +
  geom_text(aes(label = round(Importance, 2)), vjust = -0.5, size = 5) +
  theme_ipsum() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

# Guessing Score Plot
guess_plot <- ggplot(guess_data, aes(x = Variable, y = Importance, fill = Variable)) +
  geom_bar(stat = "identity", show.legend = FALSE) +
  scale_fill_manual(values = c("#C06000", "#005050")) +
  labs(
    title = "Guessing score",
    y = "Relative importance",
    x = NULL
  ) +
  geom_text(aes(label = round(Importance, 2)), vjust = -0.5, size = 5) +
  theme_ipsum() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank())

# Figure A14: Variance decomposition of misperception and guessing scores
jpeg("figure_A14.jpeg", width = 9, height = 7, units = "in", res=800)
ggarrange(mis_plot, guess_plot)
dev.off()

## Figure 4: Association of guessing and misperceptions with anti-immigration attitudes
jpeg("figure4.jpeg", width = 5.5, height = 4, units = "in", res=800)
ggplot(mosaich, aes(x = immigration)) + 
  geom_smooth(aes(y = misperception, color = "Misperception"), method = "loess", se = TRUE) +  # LOESS for misperception
  geom_smooth(aes(y = guessing, color = "Guessing"), method = "loess", se = TRUE) +  # LOESS for guessing
  scale_color_manual(values = c("Misperception" = "#005050", "Guessing" = "#C06000"),name = NULL) +  # Custom colors
  labs(  
    x = "Anti-immigration attitudes",
    y = "Score",
    color = "Legend"
  ) +
  theme_ipsum()  
dev.off()


## gender analysis

# Accuracy
t_accuracy <- t.test(inaccuracy ~ gender, data = mosaich)
print(t_accuracy)

# Confidence
t_confidence <- t.test(confidence ~ gender, data = mosaich)
print(t_confidence)

t_misperception <- t.test(misperception ~ gender, data = mosaich)
print(t_misperception)

# Confidence
t_guessing <- t.test(guessing ~ gender, data = mosaich)
print(t_guessing)

# Boxplot for Accuracy
theme_set(theme_ipsum(base_size = 14))

mosaich$gender <- factor(mosaich$gender,
                         levels = c(1, 2),
                         labels = c("Male", "Female"))

# Helper function to create a ggplot boxplot with t-test
set.seed(123)  # Set seed for reproducibility
plot_ttest <- function(data, yvar, ylab, label_y = NULL) {
  ggplot(data %>% filter(gender %in% c("Male", "Female")),
         aes(x = gender, y = .data[[yvar]], fill = gender)) +
    geom_boxplot(alpha = 0.7, outlier.shape = NA) +
    geom_jitter(width = 0.2, alpha = 0.1) +
    stat_compare_means(method = "t.test", label.y = label_y, size = 4) +
    labs(title = paste(ylab, "by gender"), x = NULL, y = ylab) +
    scale_fill_manual(values = c("#C06000", "#005050")) +
    theme_ipsum(base_size = 14) +
    theme(
      legend.position = "none",
      panel.background = element_rect(fill = "white", color = NA),
      plot.background = element_rect(fill = "white", color = NA)
    )
}

# Adjust y-axis position for t-test annotation (based on data ranges)
p1 <- plot_ttest(mosaich, "inaccuracy", "Inaccuracy", label_y = 65)
p2 <- plot_ttest(mosaich, "confidence", "Confidence", label_y = 0.95)
p3 <- plot_ttest(mosaich, "misperception", "Misperception", label_y = 45)
p4 <- plot_ttest(mosaich, "guessing", "Guessing", label_y = 65)

# Create the combined plot
combined_plot <- ggpubr::ggarrange(p1, p2, p3, p4, ncol = 4, nrow = 1)

# Save as high-resolution PNG
ggsave("figure_A18.png", combined_plot,
       width = 15, height = 5, dpi = 300)

# Ensure gender is a factor
mosaich$gender <- factor(mosaich$gender)

# Drop NAs before computing Cohen's d
library(effsize)
cohen.d(confidence ~ gender, data = na.omit(mosaich[, c("confidence", "gender")]))

cohen.d(inaccuracy ~ gender, data = na.omit(mosaich[, c("inaccuracy", "gender")]))

## Figure A13: Normalized misperception and guessing by immigration attitudes
jpeg("figureA13.jpeg", width = 6.5, height = 5, units = "in", res=800)
ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = misperception_norm, color = "Misperception"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = guessing_norm, color = "Guessing"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Misperception" = "#005050", "Guessing" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    color = "Legend"
  ) +
  theme_ipsum()  
dev.off()

## correlation plots
myvars <- c("misperception1", "misperception2", "misperception3", "misperception4", "misperception5", "misperception6", "misperception7", "misperception8")
perceptions <- mosaich[myvars]

myvars <- c("guessing1", "guessing2", "guessing3", "guessing4", "guessing5", "guessing6", "guessing7", "guessing8")
guesses <- mosaich[myvars]

# integrated correlation plot
# Extract misperception and guessing score matrices
misperception_matrix <- cor(perceptions, method="spearman", use = "pairwise.complete.obs")
guessing_matrix <- cor(guesses, method="spearman", use = "pairwise.complete.obs")

# Create a combined matrix filled with NAs
combined_matrix <- matrix(NA, nrow = 8, ncol = 8)
rownames(combined_matrix) <- colnames(combined_matrix) <- paste0("Item ", 1:8)

# Fill the upper triangle with misperception correlations
combined_matrix[upper.tri(combined_matrix)] <- misperception_matrix[upper.tri(misperception_matrix)]

# Fill the lower triangle with guessing correlations
combined_matrix[lower.tri(combined_matrix)] <- guessing_matrix[lower.tri(guessing_matrix)]

# Define new item labels
new_labels <- c("Immigrants", "Refugees", "Muslims", "Europeans", 
                "Unemployed", "Educated", "Criminals", "Undocumented")
colnames(combined_matrix) <- new_labels
rownames(combined_matrix) <- new_labels

jpeg("figure_A11.jpeg", width = 9, height = 8, units = "in", res=800)
par(mar = c(10, 10, 4, 10), mgp = c(2, 7, 4), tck = -.01, las = 1, bty = "l")  
par(mfrow = c(1, 1))  
friendly_colors <- colorRampPalette(c("#005050", "#f7f7f7", "#C06000"))(200)
# Create the correlation plot
corrplot(combined_matrix, method = "color", type = "full",
         order = 'original', diag = FALSE,
         col = friendly_colors,  # Apply friendly color scheme
         tl.col = "black", tl.srt = 45,  # Rotate labels
         tl.cex = 0.8,  # Adjust text size
         addCoef.col = "black", number.cex = 0.7,  # Show correlation values
         cl.pos = "b")  # Move legend to bottom  

# Add rotated labels for "Misperception" (upper triangle) and "Guessing" (lower triangle)
text(x = -1.3, y = 1.7, labels = "Guessing", col = "black", cex = 1.5, font = 2, srt = 90)  # 
text(x = 9., y = 7, labels = "Misperception", col = "black", cex = 1.5, font = 2, srt = 270)  # 
dev.off()

myvars <- c("scmisperception1_norm", "scmisperception2_norm", "scmisperception3_norm", "scmisperception4_norm", "scmisperception5_norm", "scmisperception6_norm", "scmisperception7_norm", "scmisperception8_norm")
perceptions <- mosaich[myvars]

myvars <- c("scguess1_norm", "scguess2_norm", "scguess3_norm", "scguess4_norm", "scguess5_norm", "scguess6_norm", "scguess7_norm", "scguess8_norm")
guesses <- mosaich[myvars]

# integrated correlation plot
# Extract misperception and guessing score matrices
misperception_matrix <- cor(perceptions, method="spearman", use = "pairwise.complete.obs")
guessing_matrix <- cor(guesses, method="spearman", use = "pairwise.complete.obs")

# Create a combined matrix filled with NAs
combined_matrix <- matrix(NA, nrow = 8, ncol = 8)
rownames(combined_matrix) <- colnames(combined_matrix) <- paste0("Item ", 1:8)

# Fill the upper triangle with misperception correlations
combined_matrix[upper.tri(combined_matrix)] <- misperception_matrix[upper.tri(misperception_matrix)]

# Fill the lower triangle with guessing correlations
combined_matrix[lower.tri(combined_matrix)] <- guessing_matrix[lower.tri(guessing_matrix)]

# Define new item labels
new_labels <- c("Immigrants", "Refugees", "Muslims", "Europeans", 
                "Unemployed", "Educated", "Criminals", "Undocumented")
colnames(combined_matrix) <- new_labels
rownames(combined_matrix) <- new_labels

jpeg("figure_A12.jpeg", width = 9, height = 8, units = "in", res=800)
par(mar = c(10, 10, 4, 10), mgp = c(2, 7, 4), tck = -.01, las = 1, bty = "l")  
par(mfrow = c(1, 1))  
friendly_colors <- colorRampPalette(c("#005050", "#f7f7f7", "#C06000"))(200)
# Create the correlation plot
corrplot(combined_matrix, method = "color", type = "full",
         order = 'original', diag = FALSE,
         col = friendly_colors,  # 
         tl.col = "black", tl.srt = 45,  # 
         tl.cex = 0.8,  # 
         addCoef.col = "black", number.cex = 0.7,  # 
         cl.pos = "b")  #
text(x = -1.3, y = 1.7, labels = "Guessing", col = "black", cex = 1.5, font = 2, srt = 90)  # 
text(x = 9., y = 7, labels = "Misperception", col = "black", cex = 1.5, font = 2, srt = 270)  # 
dev.off()

## Time stamp analysis
cor(mosaich$t_NICS8a_wp2_Page_Submit, mosaich$confidence1, use="complete.obs") #
cor(mosaich$t_NICS9a_wp2_Page_Submit, mosaich$confidence2, use="complete.obs") # 
cor(mosaich$t_NICS10a_wp2_Page_Submit, mosaich$confidence3, use="complete.obs") # 
cor(mosaich$t_NICS11a_wp2_Page_Submit, mosaich$confidence4, use="complete.obs") #
cor(mosaich$t_NICS12a_wp2_Page_Submit, mosaich$confidence5, use="complete.obs") #
cor(mosaich$t_NICS13a_wp2_Page_Submit, mosaich$confidence6, use="complete.obs") # 
cor(mosaich$t_NICS14a_wp2_Page_Submit, mosaich$confidence7, use="complete.obs") # 
cor(mosaich$t_NICS15a_wp2_Page_Submit, mosaich$confidence8, use="complete.obs") #
mosaich$response_time <- mosaich$t_NICS8a_wp2_Page_Submit + mosaich$t_NICS9a_wp2_Page_Submit + mosaich$t_NICS10a_wp2_Page_Submit + mosaich$t_NICS11a_wp2_Page_Submit + mosaich$t_NICS12a_wp2_Page_Submit + mosaich$t_NICS13a_wp2_Page_Submit + mosaich$t_NICS14a_wp2_Page_Submit + mosaich$t_NICS15a_wp2_Page_Submit
cor(mosaich$response_time, mosaich$confidence, use="complete.obs") # 

# Define custom colors
point_color <- "#005050"
loess_color <- "#C06000"

# Create the plots with jittered points and LOESS curve
set.seed(123)  # Set seed for reproducibility
plot_rt1 <- ggplot(mosaich, aes(x = log(t_NICS8a_wp2_Page_Submit), y = confidence1)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Immigrants")

set.seed(123)  # Set seed for reproducibility
plot_rt2 <- ggplot(mosaich, aes(x = log(t_NICS9a_wp2_Page_Submit), y = confidence2)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Refugees")

set.seed(123)  # Set seed for reproducibility
plot_rt3 <- ggplot(mosaich, aes(x = log(t_NICS10a_wp2_Page_Submit), y = confidence3)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Muslims")

set.seed(123)  # Set seed for reproducibility
plot_rt4 <- ggplot(mosaich, aes(x = log(t_NICS11a_wp2_Page_Submit), y = confidence4)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Europeans")

set.seed(123)  # Set seed for reproducibility
plot_rt5 <- ggplot(mosaich, aes(x = log(t_NICS12a_wp2_Page_Submit), y = confidence5)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Unemployed")

set.seed(123)  # Set seed for reproducibility
plot_rt6 <- ggplot(mosaich, aes(x = log(t_NICS13a_wp2_Page_Submit), y = confidence6)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Educated")

set.seed(123)  # Set seed for reproducibility
plot_rt7 <- ggplot(mosaich, aes(x = log(t_NICS14a_wp2_Page_Submit), y = confidence7)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Criminals")

set.seed(123)  # Set seed for reproducibility
plot_rt8 <- ggplot(mosaich, aes(x = log(t_NICS15a_wp2_Page_Submit), y = confidence8)) +
  theme_ipsum() +
  geom_jitter(color = point_color, size = 3, alpha = 0.3, width = 0.1, height = 0.1) +
  stat_smooth(method = "loess", fullrange = TRUE, color = loess_color) +
  xlab("Response time (log)") +
  ylab("Confidence") +
  ggtitle("Undocumented")

jpeg("figure_A17.jpeg", width = 12, height = 6, units = "in", res=800)
ggarrange(plot_rt1,plot_rt2,plot_rt3,plot_rt4,plot_rt5,plot_rt6,plot_rt7,plot_rt8,ncol=4,nrow=2)
dev.off()

response_time1 <- lm(log(response_time) ~ inaccuracy + confidence, data = mosaich)
summary(response_time1) # 
response_time2 <- lm(log(response_time) ~ inaccuracy*confidence, data = mosaich)
summary(response_time2) # 
response_time3 <- lm(misperception ~ inaccuracy*log(response_time), data = mosaich)
summary(response_time3) # 

mosaich_rt <- subset(mosaich, mosaich$response_time<1000)

# create alternative misperception-score
mosaich$misperception_rt <- mosaich$inaccuracy*mosaich$response_time
mosaich$guessing_rt <- mosaich$inaccuracy*((mosaich$response_time*(-1))+3802.05)

## Replication of Deutschmann (2024)

### implement method by Deutschmann (2024)

mosaich$perc1round <- mosaich$NICS8a_wp2
mosaich$perc1round[mosaich$perc1round!=0  & mosaich$perc1round!=5 & mosaich$perc1round!=10  & mosaich$perc1round!=15 & mosaich$perc1round!=20  & mosaich$perc1round!=25 & mosaich$perc1round!=30  & mosaich$perc1round!=35 & mosaich$perc1round!=40  & mosaich$perc1round!=45 & mosaich$perc1round!=50  & mosaich$perc1round!=55 & mosaich$perc1round!=60  & mosaich$perc1round!=65 & mosaich$perc1round!=70  & mosaich$perc1round!=75 & mosaich$perc1round!=80  & mosaich$perc1round!=85 & mosaich$perc1round!=90  & mosaich$perc1round!=95 & mosaich$perc1round!=100] <- NA
mosaich$perc1notround <- mosaich$NICS8a_wp2
mosaich$perc1notround[mosaich$perc1notround==0  | mosaich$perc1notround==5 | mosaich$perc1notround==10  | mosaich$perc1notround==15 | mosaich$perc1notround==20  | mosaich$perc1notround==25 | mosaich$perc1notround==30  | mosaich$perc1notround==35 | mosaich$perc1notround==40  | mosaich$perc1notround==45 | mosaich$perc1notround==50  | mosaich$perc1notround==55 | mosaich$perc1notround==60  | mosaich$perc1notround==65 | mosaich$perc1notround==70  | mosaich$perc1notround==75 | mosaich$perc1notround==80  | mosaich$perc1notround==85 | mosaich$perc1notround==90  | mosaich$perc1notround==95 | mosaich$perc1notround==100] <- NA
mosaich$perc2round <- mosaich$NICS9a_wp2
mosaich$perc2round[mosaich$perc2round!=0  & mosaich$perc2round!=5 & mosaich$perc2round!=10  & mosaich$perc2round!=15 & mosaich$perc2round!=20  & mosaich$perc2round!=25 & mosaich$perc2round!=30  & mosaich$perc2round!=35 & mosaich$perc2round!=40  & mosaich$perc2round!=45 & mosaich$perc2round!=50  & mosaich$perc2round!=55 & mosaich$perc2round!=60  & mosaich$perc2round!=65 & mosaich$perc2round!=70  & mosaich$perc2round!=75 & mosaich$perc2round!=80  & mosaich$perc2round!=85 & mosaich$perc2round!=90  & mosaich$perc2round!=95 & mosaich$perc2round!=100] <- NA
mosaich$perc2notround <- mosaich$NICS9a_wp2
mosaich$perc2notround[mosaich$perc2notround==0  | mosaich$perc2notround==5 | mosaich$perc2notround==10  | mosaich$perc2notround==15 | mosaich$perc2notround==20  | mosaich$perc2notround==25 | mosaich$perc2notround==30  | mosaich$perc2notround==35 | mosaich$perc2notround==40  | mosaich$perc2notround==45 | mosaich$perc2notround==50  | mosaich$perc2notround==55 | mosaich$perc2notround==60  | mosaich$perc2notround==65 | mosaich$perc2notround==70  | mosaich$perc2notround==75 | mosaich$perc2notround==80  | mosaich$perc2notround==85 | mosaich$perc2notround==90  | mosaich$perc2notround==95 | mosaich$perc2notround==100] <- NA
mosaich$perc3round <- mosaich$NICS10a_wp2
mosaich$perc3round[mosaich$perc3round!=0  & mosaich$perc3round!=5 & mosaich$perc3round!=10  & mosaich$perc3round!=15 & mosaich$perc3round!=20  & mosaich$perc3round!=25 & mosaich$perc3round!=30  & mosaich$perc3round!=35 & mosaich$perc3round!=40  & mosaich$perc3round!=45 & mosaich$perc3round!=50  & mosaich$perc3round!=55 & mosaich$perc3round!=60  & mosaich$perc3round!=65 & mosaich$perc3round!=70  & mosaich$perc3round!=75 & mosaich$perc3round!=80  & mosaich$perc3round!=85 & mosaich$perc3round!=90  & mosaich$perc3round!=95 & mosaich$perc3round!=100] <- NA
mosaich$perc3notround <- mosaich$NICS10a_wp2
mosaich$perc3notround[mosaich$perc3notround==0  | mosaich$perc3notround==5 | mosaich$perc3notround==10  | mosaich$perc3notround==15 | mosaich$perc3notround==20  | mosaich$perc3notround==25 | mosaich$perc3notround==30  | mosaich$perc3notround==35 | mosaich$perc3notround==40  | mosaich$perc3notround==45 | mosaich$perc3notround==50  | mosaich$perc3notround==55 | mosaich$perc3notround==60  | mosaich$perc3notround==65 | mosaich$perc3notround==70  | mosaich$perc3notround==75 | mosaich$perc3notround==80  | mosaich$perc3notround==85 | mosaich$perc3notround==90  | mosaich$perc3notround==95 | mosaich$perc3notround==100] <- NA
mosaich$perc4round <- mosaich$NICS11a_wp2
mosaich$perc4round[mosaich$perc4round!=0  & mosaich$perc4round!=5 & mosaich$perc4round!=10  & mosaich$perc4round!=15 & mosaich$perc4round!=20  & mosaich$perc4round!=25 & mosaich$perc4round!=30  & mosaich$perc4round!=35 & mosaich$perc4round!=40  & mosaich$perc4round!=45 & mosaich$perc4round!=50  & mosaich$perc4round!=55 & mosaich$perc4round!=60  & mosaich$perc4round!=65 & mosaich$perc4round!=70  & mosaich$perc4round!=75 & mosaich$perc4round!=80  & mosaich$perc4round!=85 & mosaich$perc4round!=90  & mosaich$perc4round!=95 & mosaich$perc4round!=100] <- NA
mosaich$perc4notround <- mosaich$NICS11a_wp2
mosaich$perc4notround[mosaich$perc4notround==0  | mosaich$perc4notround==5 | mosaich$perc4notround==10  | mosaich$perc4notround==15 | mosaich$perc4notround==20  | mosaich$perc4notround==25 | mosaich$perc4notround==30  | mosaich$perc4notround==35 | mosaich$perc4notround==40  | mosaich$perc4notround==45 | mosaich$perc4notround==50  | mosaich$perc4notround==55 | mosaich$perc4notround==60  | mosaich$perc4notround==65 | mosaich$perc4notround==70  | mosaich$perc4notround==75 | mosaich$perc4notround==80  | mosaich$perc4notround==85 | mosaich$perc4notround==90  | mosaich$perc4notround==95 | mosaich$perc4notround==100] <- NA

mosaich$perc1round <- mosaich$NICS8a_wp2
mosaich$perc1round[mosaich$perc1round!=0  & mosaich$perc1round!=10  & mosaich$perc1round!=20  & mosaich$perc1round!=30  & mosaich$perc1round!=40  & mosaich$perc1round!=50  & mosaich$perc1round!=60  & mosaich$perc1round!=70  & mosaich$perc1round!=80  & mosaich$perc1round!=90  & mosaich$perc1round!=100] <- NA
mosaich$perc1notround <- mosaich$NICS8a_wp2
mosaich$perc1notround[mosaich$perc1notround==0  | mosaich$perc1notround==10  | mosaich$perc1notround==20  | mosaich$perc1notround==30  | mosaich$perc1notround==40  | mosaich$perc1notround==50  | mosaich$perc1notround==60  | mosaich$perc1notround==70  | mosaich$perc1notround==80  | mosaich$perc1notround==90  | mosaich$perc1notround==100] <- NA
mosaich$perc2round <- mosaich$NICS9a_wp2
mosaich$perc2round[mosaich$perc2round!=0  & mosaich$perc2round!=10  & mosaich$perc2round!=20  & mosaich$perc2round!=30  & mosaich$perc2round!=40  & mosaich$perc2round!=50  & mosaich$perc2round!=60  & mosaich$perc2round!=70  & mosaich$perc2round!=80  & mosaich$perc2round!=90  & mosaich$perc2round!=100] <- NA
mosaich$perc2notround <- mosaich$NICS9a_wp2
mosaich$perc2notround[mosaich$perc2notround==0  | mosaich$perc2notround==10  | mosaich$perc2notround==20  | mosaich$perc2notround==30  | mosaich$perc2notround==40  | mosaich$perc2notround==50  | mosaich$perc2notround==60  | mosaich$perc2notround==70  | mosaich$perc2notround==80  | mosaich$perc2notround==90  | mosaich$perc2notround==100] <- NA
mosaich$perc3round <- mosaich$NICS10a_wp2
mosaich$perc3round[mosaich$perc3round!=0  & mosaich$perc3round!=10  & mosaich$perc3round!=20  & mosaich$perc3round!=30  & mosaich$perc3round!=40  & mosaich$perc3round!=50  & mosaich$perc3round!=60  & mosaich$perc3round!=70  & mosaich$perc3round!=80  & mosaich$perc3round!=90  & mosaich$perc3round!=100] <- NA
mosaich$perc3notround <- mosaich$NICS10a_wp2
mosaich$perc3notround[mosaich$perc3notround==0  | mosaich$perc3notround==10  | mosaich$perc3notround==20  | mosaich$perc3notround==30  | mosaich$perc3notround==40  | mosaich$perc3notround==50  | mosaich$perc3notround==60  | mosaich$perc3notround==70  | mosaich$perc3notround==80  | mosaich$perc3notround==90  | mosaich$perc3notround==100] <- NA
mosaich$perc4round <- mosaich$NICS11a_wp2
mosaich$perc4round[mosaich$perc4round!=0  & mosaich$perc4round!=10  & mosaich$perc4round!=20  & mosaich$perc4round!=30  & mosaich$perc4round!=40  & mosaich$perc4round!=50  & mosaich$perc4round!=60  & mosaich$perc4round!=70  & mosaich$perc4round!=80  & mosaich$perc4round!=90  & mosaich$perc4round!=100] <- NA
mosaich$perc4notround <- mosaich$NICS11a_wp2
mosaich$perc4notround[mosaich$perc4notround==0  | mosaich$perc4notround==10  | mosaich$perc4notround==20  | mosaich$perc4notround==30  | mosaich$perc4notround==40  | mosaich$perc4notround==50  | mosaich$perc4notround==60  | mosaich$perc4notround==70  | mosaich$perc4notround==80  | mosaich$perc4notround==90  | mosaich$perc4notround==100] <- NA

mosaich$perc5round <- mosaich$NICS12a_wp2
# mosaich$perc5round[mosaich$perc5round!="0" & mosaich$perc5round!="10" & mosaich$perc5round!="20" & mosaich$perc5round!="30" & mosaich$perc5round!="40" & mosaich$perc5round!="50" & mosaich$perc5round!="60" & mosaich$perc5round!="70" & mosaich$perc5round!="80" & mosaich$perc5round!="90" & mosaich$perc5round!="100"] <- NA
mosaich$perc5round[mosaich$perc5round!="0"  & mosaich$perc5round!="5" & mosaich$perc5round!="10"  & mosaich$perc5round!="15" & mosaich$perc5round!="20"  & mosaich$perc5round!="25" & mosaich$perc5round!="30"  & mosaich$perc5round!="35" & mosaich$perc5round!="40"  & mosaich$perc5round!="45" & mosaich$perc5round!="50"  & mosaich$perc5round!="55" & mosaich$perc5round!="60"  & mosaich$perc5round!="65" & mosaich$perc5round!="70"  & mosaich$perc5round!="75" & mosaich$perc5round!="80"  & mosaich$perc5round!="85" & mosaich$perc5round!="90"  & mosaich$perc5round!="95" & mosaich$perc5round!="100"] <- NA
mosaich$perc5round <- as.numeric(mosaich$perc5round)
mosaich$perc5notround <- mosaich$NICS12a_wp2
# mosaich$perc5notround[mosaich$perc5notround=="0" | mosaich$perc5notround=="10" | mosaich$perc5notround=="20" | mosaich$perc5notround=="30" | mosaich$perc5notround=="40" | mosaich$perc5notround=="50" | mosaich$perc5notround=="60" | mosaich$perc5notround=="70" | mosaich$perc5notround=="80" | mosaich$perc5notround=="90" | mosaich$perc5notround=="100"] <- NA
mosaich$perc5notround[mosaich$perc5notround=="0"  | mosaich$perc5notround=="5" | mosaich$perc5notround=="10"  | mosaich$perc5notround=="15" | mosaich$perc5notround=="20"  | mosaich$perc5notround=="25" | mosaich$perc5notround=="30"  | mosaich$perc5notround=="35" | mosaich$perc5notround=="40"  | mosaich$perc5notround=="45" | mosaich$perc5notround=="50"  | mosaich$perc5notround=="55" | mosaich$perc5notround=="60"  | mosaich$perc5notround=="65" | mosaich$perc5notround=="70"  | mosaich$perc5notround=="75" | mosaich$perc5notround=="80"  | mosaich$perc5notround=="85" | mosaich$perc5notround=="90"  | mosaich$perc5notround=="95" | mosaich$perc5notround=="100"] <- NA
mosaich$perc5notround <- as.numeric(mosaich$perc5notround)

mosaich$perc6round <- mosaich$NICS13a_wp2
# mosaich$perc6round[mosaich$perc6round!="0" & mosaich$perc6round!="10" & mosaich$perc6round!="20" & mosaich$perc6round!="30" & mosaich$perc6round!="40" & mosaich$perc6round!="50" & mosaich$perc6round!="60" & mosaich$perc6round!="70" & mosaich$perc6round!="80" & mosaich$perc6round!="90" & mosaich$perc6round!="100"] <- NA
mosaich$perc6round[mosaich$perc6round!="0"  & mosaich$perc6round!="5" & mosaich$perc6round!="10"  & mosaich$perc6round!="15" & mosaich$perc6round!="20"  & mosaich$perc6round!="25" & mosaich$perc6round!="30"  & mosaich$perc6round!="35" & mosaich$perc6round!="40"  & mosaich$perc6round!="45" & mosaich$perc6round!="50"  & mosaich$perc6round!="55" & mosaich$perc6round!="60"  & mosaich$perc6round!="65" & mosaich$perc6round!="70"  & mosaich$perc6round!="75" & mosaich$perc6round!="80"  & mosaich$perc6round!="85" & mosaich$perc6round!="90"  & mosaich$perc6round!="95" & mosaich$perc6round!="100"] <- NA
mosaich$perc6round <- as.numeric(mosaich$perc6round)
mosaich$perc6notround <- mosaich$NICS13a_wp2
# mosaich$perc6notround[mosaich$perc6notround=="0" | mosaich$perc6notround=="10" | mosaich$perc6notround=="20" | mosaich$perc6notround=="30" | mosaich$perc6notround=="40" | mosaich$perc6notround=="50" | mosaich$perc6notround=="60" | mosaich$perc6notround=="70" | mosaich$perc6notround=="80" | mosaich$perc6notround=="90" | mosaich$perc6notround=="100"] <- NA
mosaich$perc6notround[mosaich$perc6notround=="0"  | mosaich$perc6notround=="5" | mosaich$perc6notround=="10"  | mosaich$perc6notround=="15" | mosaich$perc6notround=="20"  | mosaich$perc6notround=="25" | mosaich$perc6notround=="30"  | mosaich$perc6notround=="35" | mosaich$perc6notround=="40"  | mosaich$perc6notround=="45" | mosaich$perc6notround=="50"  | mosaich$perc6notround=="55" | mosaich$perc6notround=="60"  | mosaich$perc6notround=="65" | mosaich$perc6notround=="70"  | mosaich$perc6notround=="75" | mosaich$perc6notround=="80"  | mosaich$perc6notround=="85" | mosaich$perc6notround=="90"  | mosaich$perc6notround=="95" | mosaich$perc6notround=="100"] <- NA
mosaich$perc6notround <- as.numeric(mosaich$perc6notround)

mosaich$perc7round <- mosaich$NICS14a_wp2
# mosaich$perc7round[mosaich$perc7round!="0" & mosaich$perc7round!="10" & mosaich$perc7round!="20" & mosaich$perc7round!="30" & mosaich$perc7round!="40" & mosaich$perc7round!="50" & mosaich$perc7round!="60" & mosaich$perc7round!="70" & mosaich$perc7round!="80" & mosaich$perc7round!="90" & mosaich$perc7round!="100"] <- NA
mosaich$perc7round[mosaich$perc7round!="0"  & mosaich$perc7round!="5" & mosaich$perc7round!="10"  & mosaich$perc7round!="15" & mosaich$perc7round!="20"  & mosaich$perc7round!="25" & mosaich$perc7round!="30"  & mosaich$perc7round!="35" & mosaich$perc7round!="40"  & mosaich$perc7round!="45" & mosaich$perc7round!="50"  & mosaich$perc7round!="55" & mosaich$perc7round!="60"  & mosaich$perc7round!="65" & mosaich$perc7round!="70"  & mosaich$perc7round!="75" & mosaich$perc7round!="80"  & mosaich$perc7round!="85" & mosaich$perc7round!="90"  & mosaich$perc7round!="95" & mosaich$perc7round!="100"] <- NA
mosaich$perc7round <- as.numeric(mosaich$perc7round)
mosaich$perc7notround <- mosaich$NICS14a_wp2
# mosaich$perc7notround[mosaich$perc7notround=="0" | mosaich$perc7notround=="10" | mosaich$perc7notround=="20" | mosaich$perc7notround=="30" | mosaich$perc7notround=="40" | mosaich$perc7notround=="50" | mosaich$perc7notround=="60" | mosaich$perc7notround=="70" | mosaich$perc7notround=="80" | mosaich$perc7notround=="90" | mosaich$perc7notround=="100"] <- NA
mosaich$perc7notround[mosaich$perc7notround=="0"  | mosaich$perc7notround=="5" | mosaich$perc7notround=="10"  | mosaich$perc7notround=="15" | mosaich$perc7notround=="20"  | mosaich$perc7notround=="25" | mosaich$perc7notround=="30"  | mosaich$perc7notround=="35" | mosaich$perc7notround=="40"  | mosaich$perc7notround=="45" | mosaich$perc7notround=="50"  | mosaich$perc7notround=="55" | mosaich$perc7notround=="60"  | mosaich$perc7notround=="65" | mosaich$perc7notround=="70"  | mosaich$perc7notround=="75" | mosaich$perc7notround=="80"  | mosaich$perc7notround=="85" | mosaich$perc7notround=="90"  | mosaich$perc7notround=="95" | mosaich$perc7notround=="100"] <- NA
mosaich$perc7notround <- as.numeric(mosaich$perc7notround)

mosaich$perc8round <- mosaich$NICS15a_wp2
# mosaich$perc8round[mosaich$perc8round!="0" & mosaich$perc8round!="10" & mosaich$perc8round!="20" & mosaich$perc8round!="30" & mosaich$perc8round!="40" & mosaich$perc8round!="50" & mosaich$perc8round!="60" & mosaich$perc8round!="70" & mosaich$perc8round!="80" & mosaich$perc8round!="90" & mosaich$perc8round!="100"] <- NA
mosaich$perc8round[mosaich$perc8round!="0"  & mosaich$perc8round!="5" & mosaich$perc8round!="10"  & mosaich$perc8round!="15" & mosaich$perc8round!="20"  & mosaich$perc8round!="25" & mosaich$perc8round!="30"  & mosaich$perc8round!="35" & mosaich$perc8round!="40"  & mosaich$perc8round!="45" & mosaich$perc8round!="50"  & mosaich$perc8round!="55" & mosaich$perc8round!="60"  & mosaich$perc8round!="65" & mosaich$perc8round!="70"  & mosaich$perc8round!="75" & mosaich$perc8round!="80"  & mosaich$perc8round!="85" & mosaich$perc8round!="90"  & mosaich$perc8round!="95" & mosaich$perc8round!="100"] <- NA
mosaich$perc8round <- as.numeric(mosaich$perc8round)
mosaich$perc8notround <- mosaich$NICS15a_wp2
# mosaich$perc8notround[mosaich$perc8notround=="0" | mosaich$perc8notround=="10" | mosaich$perc8notround=="20" | mosaich$perc8notround=="30" | mosaich$perc8notround=="40" | mosaich$perc8notround=="50" | mosaich$perc8notround=="60" | mosaich$perc8notround=="70" | mosaich$perc8notround=="80" | mosaich$perc8notround=="90" | mosaich$perc8notround=="100"] <- NA
mosaich$perc8notround[mosaich$perc8notround=="0"  | mosaich$perc8notround=="5" | mosaich$perc8notround=="10"  | mosaich$perc8notround=="15" | mosaich$perc8notround=="20"  | mosaich$perc8notround=="25" | mosaich$perc8notround=="30"  | mosaich$perc8notround=="35" | mosaich$perc8notround=="40"  | mosaich$perc8notround=="45" | mosaich$perc8notround=="50"  | mosaich$perc8notround=="55" | mosaich$perc8notround=="60"  | mosaich$perc8notround=="65" | mosaich$perc8notround=="70"  | mosaich$perc8notround=="75" | mosaich$perc8notround=="80"  | mosaich$perc8notround=="85" | mosaich$perc8notround=="90"  | mosaich$perc8notround=="95" | mosaich$perc8notround=="100"] <- NA
mosaich$perc8notround <- as.numeric(mosaich$perc8notround)

mosaich$accuracy1round <- abs(mosaich$perc1round-31)
mosaich$accuracy1notround <- abs(mosaich$perc1notround-31)
mosaich$accuracy2round <- abs(mosaich$perc2round-5)
mosaich$accuracy2notround <- abs(mosaich$perc2notround-5)
mosaich$accuracy3round <- abs(mosaich$perc3round-14)
mosaich$accuracy3notround <- abs(mosaich$perc3notround-14)
mosaich$accuracy4round <- abs(mosaich$perc4round-66)
mosaich$accuracy4notround <- abs(mosaich$perc4notround-66)
mosaich$accuracy5round <- abs(mosaich$perc5round-4)
mosaich$accuracy5notround <- abs(mosaich$perc5notround-4)
mosaich$accuracy6round <- abs(mosaich$perc6round-35)
mosaich$accuracy6notround <- abs(mosaich$perc6notround-35)
mosaich$accuracy7round <- abs(mosaich$perc7round-54)
mosaich$accuracy7notround <- abs(mosaich$perc7notround-54)
mosaich$accuracy8round <- abs(mosaich$perc8round-4)
mosaich$accuracy8notround <- abs(mosaich$perc8notround-4)

mosaich$accuracy1rounddir <- mosaich$perc1round-31
mosaich$accuracy1notrounddir <- mosaich$perc1notround-31
mosaich$accuracy2rounddir <- mosaich$perc2round-5
mosaich$accuracy2notrounddir <- mosaich$perc2notround-5
mosaich$accuracy3rounddir <- mosaich$perc3round-14
mosaich$accuracy3notrounddir <- mosaich$perc3notround-14
mosaich$accuracy4rounddir <- mosaich$perc4round-66
mosaich$accuracy4notrounddir <- mosaich$perc4notround-66
mosaich$accuracy5rounddir <- mosaich$perc5round-4
mosaich$accuracy5notrounddir <- mosaich$perc5notround-4
mosaich$accuracy6rounddir <- mosaich$perc6round-35
mosaich$accuracy6notrounddir <- mosaich$perc6notround-35
mosaich$accuracy7rounddir <- mosaich$perc7round-54
mosaich$accuracy7notrounddir <- mosaich$perc7notround-54
mosaich$accuracy8rounddir <- mosaich$perc8round-4
mosaich$accuracy8notrounddir <- mosaich$perc8notround-4


## which group is more or less accurate?
mean(mosaich$accuracy1round, na.rm=T) # 11.96711
mean(mosaich$accuracy1notround, na.rm=T) # 10.72862
mean(mosaich$accuracy2round, na.rm=T) # 34.018
mean(mosaich$accuracy2notround, na.rm=T) # 13.53521
mean(mosaich$accuracy3round, na.rm=T) # 15.81818
mean(mosaich$accuracy3notround, na.rm=T) # 9.576835
mean(mosaich$accuracy4notround, na.rm=T) # 31.03006
mean(mosaich$accuracy4round, na.rm=T) # 20.47651
mean(mosaich$accuracy5notround, na.rm=T) # 5.510072
mean(mosaich$accuracy5round, na.rm=T) # 16.80947
mean(mosaich$accuracy6notround, na.rm=T) # 27.67442
mean(mosaich$accuracy6round, na.rm=T) # 17.30277
mean(mosaich$accuracy7notround, na.rm=T) # 42.15525
mean(mosaich$accuracy7round, na.rm=T) # 23.1569
mean(mosaich$accuracy8notround, na.rm=T) # 3.655271
mean(mosaich$accuracy8round, na.rm=T) # 9.934082
# not-rounders more accurate than rounders

mosaich$accuracyround <- mosaich$accuracy1round + mosaich$accuracy2round + mosaich$accuracy3round + mosaich$accuracy4round + mosaich$accuracy5round + mosaich$accuracy6round + mosaich$accuracy7round + mosaich$accuracy8round
mosaich$accuracynotround <- mosaich$accuracy1notround + mosaich$accuracy2notround + mosaich$accuracy3notround + mosaich$accuracy4notround + mosaich$accuracy5notround + mosaich$accuracy6notround + mosaich$accuracy7notround + mosaich$accuracy8notround

plot1 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy1notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # LOESS for misperception
  geom_smooth(aes(y = accuracy1round, color = "Rounders"), method = "loess", se = TRUE) +  # LOESS for guessing
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # Custom colors
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Immigrants",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot2 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy2notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy2round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  #
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Refugees",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot3 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy3notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy3round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Muslims",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot4 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy4notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy4round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Europeans",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot5 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy5notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy5round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Unemployed",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot6 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy6notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy6round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Educated",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot7 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy7notround, color = "Non-rounders"), method = "loess", se = TRUE) +  # 
  geom_smooth(aes(y = accuracy7round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Criminals",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic
plot8 <- ggplot(mosaich, aes(x = immigration)) +
  geom_smooth(aes(y = accuracy8notround, color = "Non-rounders"), method = "loess", se = TRUE) +  #
  geom_smooth(aes(y = accuracy8round, color = "Rounders"), method = "loess", se = TRUE) +  # 
  scale_color_manual(values = c("Non-rounders" = "#005050", "Rounders" = "#C06000"),name = NULL) +  # 
  labs(
    x = "Anti-immigration attitudes",
    y = "Score",
    title = "Undocumented",
    color = "Legend"
  ) +
  theme_ipsum()  # Clean aesthetic

combined_plot <- ggarrange(
  plot1, plot2, plot3, plot4,
  plot5, plot6, plot7, plot8,
  ncol = 4, nrow = 2,
  common.legend = TRUE,      
  legend = "bottom"         
)

# Save the combined plot to a file
ggsave("figure_A16.jpeg", combined_plot, width = 12, height = 6, units = "in", dpi = 800)

# compare correlations
myvars <- c("accuracy1round", "accuracy2round", "accuracy3round", "accuracy4round", "accuracy5round", "accuracy6round", "accuracy7round", "accuracy8round")
rounders <- mosaich[myvars]
myvars <- c("accuracy1notround", "accuracy2notround", "accuracy3notround", "accuracy4notround", "accuracy5notround", "accuracy6notround", "accuracy7notround", "accuracy8notround")
notrounders <- mosaich[myvars]

# integrated correlation plot
rounders_matrix <- cor(rounders, method="spearman", use = "pairwise.complete.obs")
notrounders_matrix <- cor(notrounders, method="spearman", use = "pairwise.complete.obs")

# Create a combined matrix filled with NAs
combined_matrix <- matrix(NA, nrow = 8, ncol = 8)
rownames(combined_matrix) <- colnames(combined_matrix) <- paste0("Item ", 1:8)

# Fill the upper triangle with misperception correlations
combined_matrix[upper.tri(combined_matrix)] <- notrounders_matrix[upper.tri(notrounders_matrix)]

# Fill the lower triangle with guessing correlations
combined_matrix[lower.tri(combined_matrix)] <- rounders_matrix[lower.tri(rounders_matrix)]

# Define new item labels
new_labels <- c("Immigrants", "Refugees", "Muslims", "Europeans", 
                "Unemployed", "Educated", "Criminals", "Undocumented")
colnames(combined_matrix) <- new_labels
rownames(combined_matrix) <- new_labels

jpeg("figure_A15.jpeg", width = 9, height = 8, units = "in", res=800)
par(mar = c(10, 10, 4, 10), mgp = c(2, 7, 4), tck = -.01, las = 1, bty = "l")  
par(mfrow = c(1, 1))  
friendly_colors <- colorRampPalette(c("#005050", "#f7f7f7", "#C06000"))(200)
corrplot(combined_matrix, method = "color", type = "full",
         order = 'original', diag = FALSE,
         col = friendly_colors,  # 
         tl.col = "black", tl.srt = 45,  # 
         tl.cex = 0.8,  # 
         addCoef.col = "black", number.cex = 0.7,  # 
         cl.pos = "b")  # 

# Add rotated labels for "Misperception" (upper triangle) and "Guessing" (lower triangle)
text(x = -1.3, y = 1.7, labels = "Rounders", col = "black", cex = 1.5, font = 2, srt = 90)  # 
text(x = 9., y = 7, labels = "Non-Rounders", col = "black", cex = 1.5, font = 2, srt = 270)  # 
dev.off()

# compare confidence measurements
mosaich$confdeutsch1 <- NA
mosaich$confdeutsch1[mosaich$perc1round>=0] <- 0
mosaich$confdeutsch1[mosaich$perc1notround>=0] <- 1
mosaich$confdeutsch1a <- NA
mosaich$confdeutsch1a[mosaich$perc1round1>=0] <- 0
mosaich$confdeutsch1a[mosaich$perc1notround1>=0] <- 1
mosaich$confdeutsch2 <- NA
mosaich$confdeutsch2[mosaich$perc2round>=0] <- 0
mosaich$confdeutsch2[mosaich$perc2notround>=0] <- 1
mosaich$confdeutsch3 <- NA
mosaich$confdeutsch3[mosaich$perc4round>=0] <- 0
mosaich$confdeutsch3[mosaich$perc4notround>=0] <- 1
mosaich$confdeutsch4 <- NA
mosaich$confdeutsch4[mosaich$perc4round>=0] <- 0
mosaich$confdeutsch4[mosaich$perc4notround>=0] <- 1
mosaich$confdeutsch5 <- NA
mosaich$confdeutsch5[mosaich$perc5round>=0] <- 0
mosaich$confdeutsch5[mosaich$perc5notround>=0] <- 1
mosaich$confdeutsch6 <- NA
mosaich$confdeutsch6[mosaich$perc6round>=0] <- 0
mosaich$confdeutsch6[mosaich$perc6notround>=0] <- 1
mosaich$confdeutsch7 <- NA
mosaich$confdeutsch7[mosaich$perc7round>=0] <- 0
mosaich$confdeutsch7[mosaich$perc7notround>=0] <- 1
mosaich$confdeutsch8 <- NA
mosaich$confdeutsch8[mosaich$perc8round>=0] <- 0
mosaich$confdeutsch8[mosaich$perc8notround>=0] <- 1

## latent trait analysis
library(mirt)

# Subset the accuracy and confidence variables from the mosaich dataset
accuracy_items <- mosaich[, paste0("scinaccuracy", 1:8, "_norm")]
confidence_items <- mosaich[, c("confidence1", "confidence2", "confidence3", "confidence4", "confidence5", "confidence6", "confidence7", "confidence8")]

# Running a simple IRT model with accuracy and confidence data
accuracy_model <- mirt(accuracy_items, 1)  # 
confidence_model <- mirt(confidence_items, 1)

# Summary of the latent trait model for accuracy
summary(accuracy_model)

# Summary of the latent trait model for confidence
summary(confidence_model)







